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Abstract 

Jacobi matrices are parametrized by their eigenvalues and norming con- 
stants (first coordinates of normalized eigenvectors): this coordinate system 
breaks down at reducible tridiagonal matrices. The set of real symmetric 
tridiagonal matrices with prescribed simple spectrum is a compact manifold, 
admitting an open covering by open dense sets UJ^ centered at diagonal matri- 
ces A'^, where tt spans the permutations. Bidiagonal coordinates are a variant 
of norming constants which parametrize each open set Uj{ by the Euclidean 
space. 

The reconstruction of a Jacobi matrix from inverse data is usually per- 
formed by an algorithm introduced by de Boor and Golub. In this paper we 
present a reconstruction procedure from bidiagonal coordinates and show how 
to employ it as an alternative to the de Boor-Golub algorithm. The inverse 
bidiagonal algorithm rates well in terms of speed and accuracy. 

Keywords: Jacobi matrix, inverse eigenvalue problem, bidiagonal coordinates. 
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1 Introduction 



Recall that a real tridiagonal symmetric matrix 
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is a Jacobi matrix if 6^ > for all i. Jacobi matrices have simple spectrum and 
their eigenvectors have nonzero first and last coordinates. Thus, a Jacobi matrix T 
diagonalizes uniquely as T = Q*AQ, A — diag(Ai < • • • < A„), provided we demand 
that the norming constants Wi = Qn be positive for all i. Let be the set of 
Jacobi matrices with given simple spectrum A and S"^^ — {(wi, . . . ,Wn) \ Wi > 
0,wl + - ■ • + = 1}, the open positive octant of the unit sphere S"^^ C M". Define 
the map of norming constants uj\ : S"""'^ by (jJ\{T) — w = {wi, . . . ,Wn)- 

Moser (0) proved that wa is a diffeomorphism. 

On another route, numerical analysts considered the problem of reconstructing 
a Jacobi matrix T from its eigenvalues A^, i — l...n, and the eigenvalues /i,;, 
i = 1 ... n — 1, of its bottom principal (n — 1) x (n — 1) minor. The interlacing 



*The authors acknowledge support from CNPq, IM-AGIMB and FAPERJ. 



theorem requires that Ai < /ii < A2 < /i2 < • • • < Mn-i < ^n- Existence and 
uniqueness of T, continuous dependence on A's and /x's and an iterative algorithm 
to obtain T were obtained by Hochstadt, Gray, Wilson and Hald (|^, 0, 0], 
A direct, stable algorithm was obtained by de Boor and Golub ( 1 ). The first step 
of their algorithm is the computation of the norming constants w in terms of A's 
and fi^s. The problem then boils down to computing the inverse map lijJ^^{w), as 
will be described in section 2. For a survey of the Jacobi reconstruction problem, 
see 0- 

Let be the closure of J7j, the set of tridiagonal symmetric matrices T with 
spectrum A and bi > 0. The map of norming constants wa extends continuously to 
J7a but this extension is no longer injective. Indeed, for n = 3 and A = diag(l, 2, 4) 
we have 

1 0\/lO 0\/lO 0\ 

3 — cos 2t sin 2t = cos t sin i A cos t — sin t 
sin2i 3 + cos2ty \0 — sint cost J yO sini cost J 

and therefore u}\(T) = (1,0,0), /ii ~ 4 and /i2 ~ 6 for all such T. Thus, in some 
sense, any reconstruction algorithm either from A's and /x's or from A's and w's 
must degenerate at some points of the boundary of J7a- 

In , the authors introduced hidiagonal coordinates, a variant of norming con- 
stants which behaves well at the boundary. In this paper we provide a direct recon- 
struction algorithm from bidiagonal coordinates with good behavior at boundary 
points, where norming constants break down. The conversion from norming con- 
stants to bidiagonal coordinates is simple, and the resulting algorithm is comparable 
in time and space with that of de Boor and Golub, being more accurate in many 
cases. 
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Figure 1: The manifold Xa for A = diag(l, 2, 4) 

Matters become clearer with some geometric vocabulary. Let Za 3 be the 
set of tridiagonal symmetric matrices with spectrum A. As proved in X\ is 
a compact oriented manifold of dimension n — 1. The closure J\ C I a of is 
homeomorphic to the convex polytope Va with vertices A'^ = diag(A7r(i), ■ • • , \(n)) 
where tt spans the set of permutations of {1, 2, . . . , n} (^Hli El)- Each of the 2"~^ 
possible choices of signs for the entries hi define a closed subset of 2a which is 
isomorphic to J7a: as is well known, dropping the signs of the off-diagonal entries 
of a tridiagonal symmetric matrix does not change its spectrum. Thus, Xa can 
be constructed by glueing 2"~^ copies of Vk along faces consisting of reducible 
tridiagonal matrices (i.e., matrices for which some hi is zero). In figure ^ we show 
what happens for n = 3. The polytope Vk is a hexagon and in each of its copies 
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we indicate the signs of bi and &2- Vertices are diagonal matrices and edges consist 
of reducible matrices; edges with the same label are glued. It follows that 2a is a 
bitorus for n = 3. 

Removal of the outer boundary edges in the picture yields an open dense subset 
of Xa centered at A. As we shall see, bidiagonal coordinates can be smoothly defined 
on this set, yielding an explicit diffeomorphism with R^. More generally, for each A 
and each permutation tt, we define an open dense subset U]l of Ia- The complement 
of UjI consists of matrices T for which there exist i < k < n with bk — and A^(i) 
belongs to the spectrum of the bottom principal (n — k) x (n — k) minor. As 
in the example above, is centered at A'^ and bidiagonal coordinates provide a 
diffeomorphism between and Also, the sets form an open cover of 2\: 

this is the crucial property for the local study of iterations preserving tridiagonality 
and spectrum as performed in [H]. 



2 The de Boor-Golub algorithm 



In this section we present part of the contents of 1 phrased in such a way as to 
emphasize the differences and similarities between this more well known algorithm 
and the inverse bidiagonal algorithm, to be presented in section 4. We assume that 
the off-diagonal entries bi of T are positive so that T is a Jacobi matrix. Write 
T = Q*AQ where A = diag(Ai, . . . , A„) (in this section, the simple eigenvalues Ai 
are taken in an arbitrary order), and Q is orthogonal with positive first column. 
For Db = diag(l,6i,6i62, ■ •^,&i&2 • • -^n-i) and = diag(Qii, (32i, • ■ • , Qni) = 
diag(ti;i, W2, ■ ■ ■ , Wn), write P = Di,Q*D^^ so that 
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Let pI._i = e*j.P be the fc-th row of P; in particular, 
TP = PA and the rows of P satisfy the recursion 



= PAP" 



(1,1,...,!). We have 



Pi =Po^-aiPo, 
Pk+i = P*k^ ~ o,k+ipl - blpl_j^, < A: < n - 1, 

Furthermore, the vectors pk form an orthogonal basis under the inner product 
{{u,v)) — (u, D^w). For < fc < n, let pk be the unique polynomial of degree less 
than n satisfying Pfe(Aj) — {pk)j- we have po — 1, 

Pi^tpo-aipo, pk+i ^ tpk ~ ak+iPk - blpk-i, < k < n, 

so that Pk is a monic polynomial of degree k. 

In the notation of [Jj, let Pk{t) = det{tl — Tk) where Tk is the principal minor 
of T consisting of the first k rows and columns; set also pa — 1. The expansion of 
the determinant along the last row of each minor yields 



Pi = tpo - aipo, pk+i = tpk - ak+iPk - hkPk-i 



< k < n 



(1) 



and therefore pk = pk since po = po and both sequences satisfy the same recurrence. 
Equivalently, the j-th coordinate of the vector is pk{Xj)- 
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Summing up, assume A and Dm — diag(wi, . . . ,Wn) given. The linear bijection 
between the space of real polynomials of degree less than n and M" given by eval- 
uation on the Aj's allows us to pull back the inner product ((■, •)) giving rise to an 
inner product on polynomials: 

n 

((91,92)) = ^w;|gi(Aj)(72(Aj). 

i=i 

The de Boor-Golub algorithm now constructs the monic polynomials pk using the 
orthogonality condition {(pk , Pk' )) = for k ^ k' . Recursion ^ for the polynomials 
Pk obtains the entries of T. 



3 Bidiagonal coordinates 



We quote some of the results and notations of [H] to be used in the inverse bidiagonal 
algorithm. Diagonahze T G Xa as T = Q*AQ, A = diag(Ai < A2 < • • • < A„), and 
factor Q = PLU where P is a permutation matrix, L is lower unipotent and U 
is upper triangular. For a permutation tt, let P^r be the permutation matrix with 



PttiPtts and P^e^ 



e,(,)). The PLU 



(i-ij) entry equal 1 iff i = 7r(j) (thus Pn-ir 
factorization can usually be done for several choices of the permutation matrix: it 
turns out that we can take P = P^ if and only if P G U]l, the open dense subset 
of 2a presented in the introduction. Write Qt^ = EP^^Q — L-^Ut^ where P' is a 
diagonal matrix with diagonal entries equal to 1 and —1, L^r is lower unipotent 
and C/jr is upper triangular with positive diagonal. The rows of Q^^ are eigenvectors 
of T but their first coordinates are not necessarily nonnegative: instead, signs are 
determined from the fact that the determinants of leading principal minors of 
are positive. 

Let P„ = L^^KLt, — R^^TRt, where P^r = C/^^ so that = Q-kR-k- From the 
first formula, P^ is lower triangular; from the second, it is upper Hessenberg; thus, 
P^ is lower bidiagonal: 



Ptt — 
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The map ip-^ '■ ^ M"^-'^ taking T to the n-bidiagonal coordinates {Pi , . . . , Pn-i) 
is a diffeomorphism. Indeed, start from an explicit formula for the matrix L^r in 
terms of bidiagonal coordinates: 



P, 



PI 



(AJ-AJ)(AJ-AJ) 



01 02 -01-1 



02-01-1 



\(K->-l)(K->-l)---(K-K-i) {K->-l)---(K~K-i) 







■•■ 
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Given Pjr, its QR factorization yields Q-n and P^r, from which one obtains T = 
Ptt Ptt P^^- A Straightforward computation shows that hi and /3f have the same 
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sign and that near a diagonal matrix, /3f equals bi to first order; more, for the 
inverse map <^7r = '■ K"""^ ^ ^1 C 2a, 



6^(0, . . . , 0) + D<p^{0, . . . , 0)(ui, . . . , 



Ui M2 
M2 AJ 
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For any permutation tt and any T S Z//^) the norming constants wf = WT^^i) and 
the TT-bidiagonal coordinates (3^ are related by 

• • • /?fLi 9 ^ ■ ^ 

'^(Ar-A7)---(A7-A7_,)' 



Hi 



(Aj+i - Ai ) • • • (Af^i - Af 
(A7-Ar)---(A7-A7_iX 



1< i < n - 1. 



The inverse bidiagonal algorithm, presented in the next section, obtains the 
matrix iJ^ in another way, closer in spirit to the recursions in section 2. The 
basic version of this algorithm receives as input a permutation tt, eigenvalues Af 
and bidiagonal coordinates /3f and returns the corresponding tridiagonal matrix 
T e The de Boor-Golub algorithm, instead, receives as input the eigenvalues 
A, and the norming constants Wf. in this case a simultaneous permutation tt is 
innocuous, at least with exact arithmetic. 



4 The inverse bidiagonal algorithm 

We first describe a preliminary version of the algorithm, which only works in 
the irreducible case, where all /3J (or, cquivalently, all bk) are nonzero. Write 
T = RB^R-^ where R = DbR is an upper triangular matrix with rows f^. 
Clearly, rl = elR = elDbR^ = e^R^ = elQ*L„ = (Q^ei)*L^ = {L„U„ei)*L„ = 
^11 (-t'7rGi)*-^7r — ^^'iic*-Zv^i/7T and therefore vi — cL^L-j^ ei, the value of c = un > 
being irrelevant throughout the algorithm. 
Equate rows in TR = RB^, 

\ 

A3 

to obtain f^_,_^ = f^B^ — afcf^ — b1_j^fl_^. Since R is known to be upper triangular, 
this recursion, together with the initial term f*, allows us to compute the coefficients 
ai, i — 1, . . . ,n and bf, i — l,...,n — 1. The numbers bi and Pf have the same sign 
in Uj[: this completes the preliminary version of the reconstruction algorithm for 
irreducible matrices. 

We need to modify the algorithm in order to extend it to the general case. For 
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an integer A: > 0, set 



1 



(AJ-Ar)(AJ-AJ) 



and B-^ k = l-^'^-^Tr.fe, or, more explicitly, 
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Still in the irreducible case, define = diag(l, /3f , /3J^/3J, . . . , /3J • • • and 
^ = c^^RD^^ with rows so that B^,2 = OpB-^D^^ and = ^S7r,2- Straight- 
forward computations verify that fi — L* ji^^oci- Expanding the matrix products 



as above we obtain the recursion 



fc+i 



hi- 



1' fc-i- 



Thus, from 



f^._j^ and we compute r^i?7r,2, then and and finally f'^j^^- This completes 
the description of the inverse bidiagonal algorithm for irreducible matrices; we now 
prove that this procedure works for any Z?'^ e M"~^, obtaining all matrices T G Uj[. 
Let Upper^(]R, n) be the group of upper triangular matrices with positive diagonal. 



Proposition 4.1 There is a smooth function p : lAJ^ Upper"'' ( 
p{T) = (i?7r)iiZ?t,i?7r^^^ for all irreducible matrices T G Z^J. 



I, n) satisfying 



Here, as in section 3, T = Q^K^Q-,^, L-„ = Q-^Rtt, Qtt orthogonal, L^r lower 
unipotent and Rj^ G Upper"'' (M, n) . The purpose of this proposition is to make sense 
of R for reducible matrices T (or, equivalently, for (3'^ with some zero coordinate). 
The formula in the statement defines p{T) as R for irreducible T but otherwise 
involves divisions by zero. 

Proof: Define p : Uj{ M"'^" row by row: let p^. denote the fc-th row of p{T) 
and set pi = L* jiTr.oei, pl^i — p*kB-K,2 - akPl — bl^iPl-i- The function p is 
clearly smooth in U^. Also, we proved above that p{T) = p{T) = R for irreducible 
T . Thus, by continuity, p{T) is always upper triangular with nonnegative diagonal 
entries pk,k- In the irreducible case. 



Pk,k 



(-Rtt)*;, fe- 



lt remains to prove that p^, ^ ^ for reducible T so that we can then set p — p. 

One way of completing the proof is recalling from 8 that the quotients bj/(3J 
are smooth positive functions on U^. Alternatively, from B„ 2 = L^\^L„ 2j we can 
write {B* 2)''^^pi = L^ 2^''~^L.^fiei. The coordinates of L^^qBi are all nonzero by 
equation[21and, from the standard Vandermonde argument, the vectors A'^~^ Lj^fiCi 
k = 1, . . . , n, form a basis; since L* 2 is invertible, so do the vectors {B*^2Y 
From the recursion formula, so do the vectors pk and we are done. 



ill 
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In general, we start from and use the recursive formula rl^i — r^^B^^^ ~ 
«fc^fc — bf._j^f^_^. More precisely, assume by induction that and are known. 
From the proposition, rk,k and fk-i,k-i are positive. The first nonzero coordinate 
of r^i?7r,2 occupies position fc — 1 and equals {P^_i)'^fk,k- The algorithm then 
calculates bk-i = P^-iV ^k,k / fk-i,k-i, so that the square root is evaluated at a 
strictly positive number. Notice that the algorithm treats uniformly all £ R"~^, 
i.e., there is no checking of signs or division into cases. The values of ak and of 
fl_^_i, the (fc + l)-th row of the upper triangular matrix R, are now easily obtained, 
concluding the computation of , • • ■ , Pn-i) the description of the inverse 

bidiagonal algorithm. 



5 Accuracy and tight permutations 

Empirical evidence indicates that a good choice of the permutation tt is extremely 
important for the accuracy of the inverse bidiagonal algorithm. One is reminded 
of Gaussian elimination, where pivoting strategies have a similar effect. There is 
a crucial difference, however. In Gaussian elimination, the permutation is chosen 
along the process; the inverse bidiagonal algorithm admits no easy way to acco- 
modate a change of permutation in mid- flight. As to estimating accuracy with 
respect to the choice of permutation, our theoretical understanding is limited and 
we provide instead a simple numerical experiment ^ . We start with random inverse 
data for an 8 x 8 matrix and perform the inverse bidiagonal algorithm for each 
of the 8! permutations with 8 digits of precision. Results are then compared with 
the "correct" answer, computed with an exaggerated number of digits. Different 
permutations yield very different errors: the smallest error is 2.0 • 10^^, there are 19 
other permutations with error smaller than 3 • 10^^ and there are 8 permutations 
with error greater then 7 • 10^^. The error here is defined as 

n n— 1 

\ai - a^l + ^ \bi - h\ 

i=l i=l 

where Oi and bi are the "correct" values and hi and hi are the computed values. 
All entries of T have absolute value smaller than 1. In this section, we present a 
strategy for choosing tt. 

Let Ti be the transposition (fc, fc + 1) in cycle notation. Two permutations ttq 
and TTi differ by Tk if tti = ttq ot^. Thus, K^'^ is obtained from A^" by interchanging 
the (fc, fc) and (fc -|- 1, fc -I- 1) entries. Bidiagonal coordinates (3'^° and /3f ^ are equal 
except for 



where 



<lk 



To 



N^o X^O ' 

'^k+l ~ ■^k 



as can be proved from the formulae relating /3's and w's in section 3. 

Given inverse data Xi and Wi, we call a permutation tt tight if 1 9^1 < 1 for all 
j = 1, . . . , n — 1 and we say that the transposition is tt -tightening if \q'^ \ > 1. From 
equation 131 it is easy to see that if is 7r-tightening then Ig^"^*"! < 1. Clearly, tt is 
tight if and only if there are no 7r-tightening transpositions. A tightening sequence 
is a maximal sequence (jTm) of permutations such that TTm+i differs from tt^ by a 
7r„i-tightening transposition Tk^^. Thus, a tightening sequence is either infinite or 
ends at a tight permutation. 

^ Maple worksheets for all experiments in this paper are available at 
jhttp: / /www. mat .puc-rio.br/"' nicolau/papers/invbi-mw. 
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Lemma 5.1 Tightening sequences are finite. 



Proof: Assume by contradiction that there exists an infinite tightening sequence: 
clearly, there exist mp < mi with -Kma = ""mi ■ We show that there are no such 
cycles. 

Set pk,-K — Y\i>k ^'^^ lexicographical order to define a total order in 

the permutation group: ttq -< tti if and only if there exists k' such that Pk.ixa — Pk.ni 
for k < k' and Pk' ,tto < Pk', rri- For any m, it follows from equation |21 that Pfc,7r„+i = 
Pfc,7r„ for k < km and Pfc„,,r,„+i < P/c„,7i-„, implying iTm+i -< T^m- By transitivity, 
TTmo ^ ""mi = ""mo i 8- Contradiction. ■ 

In the example discussed above, there were 4 tight permutations with errors 
between 2.8 ■ 10~^ and 5.3 • 10~^. Empirical evidence shows that this is frequent: 
tight permutations usually yield small errors in the inverse bidiagonal algorithm. 
Thus, upon receiving inverse data Xt and Wi > for a Jacobi matrix T, we first order 
the w's in decreasing order to obtain a permutation ttq and then apply tightening 
transpositions until we reach a tight permutation. Experiments suggest that this 
takes approximately n/2 sweeps. 



6 Operational costs 

We estimate the number of operations (or flops) and the amount of memory neces- 
sary to execute the inverse bidiagonal algorithm. As usual, we only keep track of 
the number of products and quotients. 

The process of finding a tight permutation will not be carefully examined: suffice 
it to say that, from empirical evidence, the number of operations is approximately 
Cn^ where C < 1/2. 

The matrices -B,r,2 and 2 will come up along the algorithm and it is therefore 
convenient to compute and keep the squares (/3f )^, with an initial cost of n opera- 
tions and n storage units. The first major step of the algorithm is the computation 
of fi = ij-^oei- For n = 4, after reordering terms, fi becomes 

\{Xl - A?)2(AJ - A-)2(AJ - AJ)2 (Af - A?)2(AJ - AJ)^ + (AJ - A?)^ + ' 

{P^)Hl3tf iP^f 1 

(AJ - AJ)(AJ - AJ)2(AJ - AJ)2 + (AJ - A7)(Af - AJ)^ + (AJ - AJ) ' 

iPil , 1 

(AJ - A?)(AJ - AJ)(AJ - AJ)2 (AJ - A^IAJ - AJ) ' 



(AJ-An(AI-AJ)(AJ-A-)y 

These terms are computed from bottom to top of each column, following the obvious 
patterns, with a cost of approximately 3n^/2 operations and n storage units. 

The recursion formula which obtains bk-i, ak and fk+i only requires fk-i and 
fk so that we only need to keep at most three rows of the triangular matrix R at 
any given time. The number of operations is approximately 2n^] also, n — 1 square 
roots are needed. 

Summing up, given a tight permutation tt and the values of /3f , a run takes 
approximately 7n^/2 operations, n square roots and 4n storage units (provided 
some units do double duty, first as entries of R and later as a's or &'s). 
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7 Benchmarks 



In this section, we compare the de Boor-Golub and inverse bidiagonal algorithms 
in a few scenarios. This is only possible for irreducible matrices since otherwise, 
as we saw, the norming constants w break down. The inverse bidiagonal algorithm 
receives as input permuted eigenvalues Af , « = 1, . . . , n, and bidiagonal coordinates 
f]^ , i = 1, . . . ,n — 1. In order to allow for comparisons, we must step back and 
provide as input the eigenvalues and the norming constants wf. we then obtain 
a tight permutation tt and compute l3f . 

It is a common feature of both algorithms that the coefficients a;, i = 1, . . . , n 
and bi, i = 1, . . . , n — I, are obtained in the order ai, bi, a2, 62, 03, • ■ ■• Also, both 
algorithms admit a reversal by conjugation. More precisely, let Pp be the permu- 
tation matrix with {Pp)ij = I if and only iii + j~n + l. Let A^ and Wi be the 
inverse data for a Jacobi matrix T: the inverse data for T = PpTPp is A.^ and 

c 

for some positive normalizing constant c (^). From data A^ and Wi either algorithm 
obtains, in this order, ai = a„, bi = bn-i, 0.2 = an-i, ■ . ■• Experiments show that it 
is far wiser to do both things, i.e., to compute the top half of T directly from Wi 
and the bottom half from Wi. In the examples below, this strategy, the two-sided 
algorithms, is always adopted. 

We implement in a Maple worksheet both the two-sided de Boor-Golub (BG) 
and the two-sided inverse bidiagonal (BI) algorithms, generate a sequence of random 
inverse problems and compare errors for different values of the dimension and of 
the number of significative digits. 

In the first class of examples, random real symmetric tridiagonal matrices T are 
obtained as follows: the nonzero entries are independent random variables with a 
Gaussian distribution centered at with variance 1. We then compute the inverse 
variables Ai and Wi of T and test the algorithms with these inputs: the norming 
constants Wi typically span several orders of magnitude. There are cases where 
either algorithm outperforms the other, but in the average BI fares decisively better 
than BG. In the worksheet, we repeated this experiment 40 times with dimension 
n = 40, working with 12 significant digits; errors were measured as in the previous 
section. The run is declared a failure if the error exceeds 0.1: there were two runs 
where both BG and BI failed, another 30 failed runs for BG and none other for BI. 

We next consider matrices near Tq, the Jacobi matrix with diagonal entries 
equal to and off-diagonal entries equal to 1: for our purposes, Tq is as good as 
the free Laplacian. It turns out that Tq is special from several points of view: the 
values of both Ai and Wi can be obtained explicitly, there are no small gaps in the 
spectrum and all norming constants Wi have roughly the same size. These features, 
particularly the last one, seem to favor BG. Indeed, for n = 40, working with 12 
digits, Gaussian perturbations of Tq with small variance in the off-diagonal entries 
favors BG: among 40 examples, there are no failures of either algorithm but the 
errors are smaller for BG than for BI. 
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